function y = outlier (x, thd, fillval)

   ## usage: y = outlier (x, thd, fillval)
   ## 
   ## detect outliers based on threshold thd, and (optionally) replace with fillval

   pkg load nan

   if nargin < 2
      thd = 100 ;
   endif
   if nargin < 3
      fillval = nan ;
   endif

   y = z = x ; I = false(size(x)) ;
   iter = ceil(thd/10) ;
   for j = 1:iter
      if sum(isfinite(x(!I))) < 10, continue ; endif
      w = x(!I) ;
      ##z(!I) = (w - min(w)) ./ (max(w) - min(w));
      z(!I) = zscore(x(!I)) ;
      I = I | abs(z) > thd ;
   endfor

   if (s = sum(I(:))) > 0
      warning("xds:xds", "outlier: found %d outliers in %d elements.", s, numel(x)) ;
      y(I) = fillval ;
   endif

endfunction
